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Abstract — We present SPH formulations of Dedner et al's 
hyperbolic/parabolic divergence cleaning scheme for magnetic 
and velocity fields. Our implementation preserves the conserva- 
tion properties of SPH which is important for stability. This is 
achieved by deriving an energy term for the ip field, and imposing 
energy conservation on the cleaning subsystem of equations. This 
necessitates use of conjugate operators for V • B and \7ip in 
the numerical equations. For both the magnetic and velocity 
fields, the average divergence error in the system is reduced by 
an order of magnitude with our cleaning algorithm. Divergence 
errors in SPMHD are maintained to < 1%, even for realistic 
3D applications with a corresponding gain in numerical stability. 
Density errors for an oscillating elliptic water drop using weakly 
compressible SPH are reduced by a factor of two. 

I. Introduction 

Magnetic fields have the property of being divergence free, 
that is V B = 0. Incompressible fluids have a similar 
divergence free property for the velocity field. Maintaining 
these divergence constraints is one of the central difficulties 
in performing accurate simulations of magnetohydrodynamics 
(MHD) and incompressible fluid behaviour. For MHD in 
particular, the presence of magnetic monopoles introduces a 
spurious force which, when large, is disruptive to the dynamics 
of the system. 

Similar approaches can be utilised to satisfy the divergence 
constraints in both cases. For example, projection methods 
construct a divergence free vector field via the solution of 
a Poisson equation and have been applied successfully to both 
systems. Specialised approaches have also been developed for 
each case. One example is the constrained transport method 
[1] for MHD, which by conserving magnetic flux through a 
closed surface, can keep the divergence constraint to within 
machine precision. For SPH simulations of incompressible 
fluids, a stiff equation of state can be used to limit density 
variations to ~ 1% [2], creating a weakly compressible fluid 
approximating incompressibility. 

The hyperbolic divergence cleaning method of Dedner et al 
{3} was introduced for maintaining the V • B = constraint in 
MHD. It involves the addition of a new scalar field, i/j, which 
is coupled to the magnetic field by 

dB 1 (i) 



dt 



This -0 field evolves according to 



and combined these produce a damped wave equation 
S 2 (V-B) ~ 2m _ , lfl(V-B) 



dt 2 



c^(V-B) 



dt 



0. 



(3) 



Thus divergence is spread away from sources by a series of 
damped waves. The wave speed, Ch, is typically chosen to be 
the fastest wave obeying the Courant stability condition. The 
damping timescale, r, acts as a diffusion on the divergence. 
By using waves to spread the divergence over a larger volume, 
the amplitude of any single large source is diminished and 
the diffusion is more effective. While originally proposed for 
use on the magnetic field for MHD simulations, this approach 
would be valid for any vector field. The damping timescale is 
set to t -1 = ach/h, where h is the smoothing length and a 
is a dimensionless quantity specifying the damping strength. 

Hyperbolic divergence cleaning has found popular use 
in both Eulerian (4), (5) and Lagrangian based codes (6), 
[7], chiefly for its simplicity, easy implementation, and low 
computational cost. However, for the SPH implementation of 
MHD (SPMHD), this method has not been widely adopted. 
Initial implementation attempts by Price (8) found divergence 
reductions were not substantial (a factor ~ 2), and the method 
risked actually increasing divergence in certain test cases. 

The work presented here describes a new formulation of 
hyperbolic divergence cleaning for SPH that removes previous 
difficulties [9]. Implementations for both the magnetic and 
velocity fields are presented. Our formulation imposes the 
constraint of energy conservation on the subsystem of cleaning 
equations, guaranteeing that energy transferred to the ip field 
must either be conserved or dissipated. This prevents increases 
in divergence. 

The paper is laid out as follows: Sec. [n] discusses hyperbolic 
divergence cleaning for the magnetic field of SPMHD. A 
brief description of SPMHD is presented (Sec. |II-A| ), along 
with the Euler Potentials (Sec. II- A 1 ) and artificial resistivity 



(Sec. II- A2) since they will be used as a basis of comparison 
for the new divergence cleaning method. In Sec. |II-B| the 



energy contained in the i/j field is derived and modifications 
are made to the cleaning equations to conserve energy, then 
the energy conserving SPMHD implementation is constructed 
(Sec. |II-C| ). Hyperbolic divergence cleaning for the velocity 
field is discussed in Sec. [TTTl Starting from an outline of 



c 2 h V-B 



(2) 



weakly compressible SPH (Sec. Ill- A ), a new energy term 
is created for the tp field for contributions from the velocity 
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field (Sec. |III-B| ) which is used to create the conservative 
SPH implementation (Sec. |III-Q . Tests of our method are 
presented in Sec. [IV| applied to three MHD problems and 



one incompressible fluid problem. The SPMHD tests include a 
simple free boundary test (Sec. |IV-A| , the Orszag-Tang vortex 
where the cleaning method is compared against resistivity and 
Euler Potentials (Sec. |IV-B| ), and a collapsing molecular cloud 
involving star formation (Sec. IV-Q . A test of the velocity 
cleaning is presented on an oscillating elliptic water drop 
using weakly compressible SPH (Sec. |IV-D| ). Conclusions are 
presented in Sec. [V] 

II. Hyperbolic divergence cleaning for the MHD 

EQUATIONS 



A. Smoothed particle magnetohydro dynamics 

The equations of ideal MHD solved in SPH are given by 

Pa = ^2rn b W ab (h a ), 

b 

l/n d im 



Pa 



(4) 
(5) 



dVq 

dt 



dB Q 

dt 



M 



VaPl 



^-V a W ab (h a ) 



^ • V h W ah (h h ) 

\ Li ^ 



^aP 



L bP b 

v ab (B a -V a W ab (h a )) 
B a (v a6 - V a W ab (h a )) 



(6) 



(7) 



Here, v is the fluid velocity, B is the magnetic field, and 
d/dt is the material derivative. The density, p, is calculated 
via summation using an iterative procedure to self consistently 
determine the smoothing length, h. Variable smoothing length 
gradients are accounted for with the ft terms (see |TQ|). 
The momentum equation contains contributions from thermal 
pressure, P, and the Lorentz force, given in terms of the 
Maxwell stress tensor 



M = BB — (P+|P 2 )l. 



(8) 



The contribution from any spurious B(V • B) force is sub- 
tracted out by including the additional term 



dt 



V B 



-B Q ^m b 



VaPl 
B b 



V a W ab (h a ) 



n 2 V a W ab (h b ) 



(9) 



The induction equation is derived from dB/dt = V x (v x B) 
with the monopole contribution removed. Hence, this scheme 
is formally equivalent to Powell's eight wave approach (TTJ. 



1) Euler Potentials: One method for maintaining the di- 
vergence constraint on the magnetic field is to use the Euler 
Potentials, defining B = Vax V/3. The potentials are advected 
exactly, representing the field lines being frozen to the fluid, 
removing the need to solve the induction equation ([7]). This has 
had reasonable success controlling divergence error to ~ 1% 
(ie, |T2| , (13)), however the Euler Potentials place limitations 
on the possible field configurations which can be represented. 
Magnetic field windings in particular cannot be modelled past 
one rotation, and such topologies would be anticipated for 
many astrophysical problems of interest. 

2) Artificial Resistivity: Artificial resistivity is added to 
SPMHD to capture magnetic shocks and discontinuities. It is 
similar to artificial viscosity, with form 

dB 



dt 



Pa ^2 mb 



Pab 



(B a -B 6 )f. V a W ab , 



(10) 



where is a dimensionless quantity of order unity, v s { g 
is a signal velocity, and p ab is the average density between 
particles a and b. This is representative of real resistivity, 

-=i|V(V-B)-i]Vx(Vx B), (11) 

and as such provides diffusion of magnetic divergence. In 
some cases, this may be sufficient to control errors, however it 
is a poor tool to control divergence error since it also dissipates 
the physical portions of the field. 

B. Hyperbolic magnetic divergence cleaning 

If just the cleaning system of equations given by ([T]) and 
^ are considered, then the total energy can be written as 

B 2 



E : 



/ 



2p p 



pdV, 



(12) 



which is the sum of magnetic energy and as yet undetermined 
energy contained in the ip field. By conservation of energy, 



pdV = 0, 



dE _ f B /dB\ de^ 
dt J pop \dtj^ dt 

and if is assumed to have differentiate form 
de^p difj 
~dt =X ldt' 

then by inserting ([T]) and ([2]), we can obtain 



/ 



B 



- x^V • B 



pdV = 0. 



(13) 



(14) 



(15) 



/ 



P0 Js 



i(jB • ds = 0. (16) 



Integrating the first term of ( [T5] ) by parts will yield 

A- X c 2 J (V-B)pdX 
pop J po 

The surface integral may be ignored. Similar terms appear, for 
example, in the SPH continuity equation which are likewise 
taken to be zero |14|. The remaining term thus implies 

X = ip/popcl and therefore 



^ 2 



2popc h 



(17) 
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Inserting this energy term into §13\ will produce 



/[-•( 

J POP \ 



dB\ 



p pc 2 h dt 2p p 2 a c 2 dt 



t/j 2 dp 

272" 



pdV = 0. 
(18) 

From the preceding analysis, it is clear that energy changes 
from the first two terms will be balanced by each other. To 
account for the third term, the evolution equation for ip can 
modified to 



difj 
~dt 



-c\V • B 



t/. 



(19) 



replacing ([2]). 



C. Discretised hyperbolic magnetic divergence cleaning 

Hyperbolic divergence cleaning is implemented into 
SPMHD using the differenced derivative operator for V • B, 



V B a = -— !— V m h (B a - B 6 ) • VW ah {h a ). 

il a Pr ^ 



(20) 



Other operator choices are permissible. It may seem that using 
the same operator as in the momentum equation (see ([9])) 
would be desirable, however we have found that doing so leads 
to excessive magnetic energy dissipation. This occurs because 
that operator also measures the disorder in the particle arrange- 
ment, which the cleaning method attempts to compensate for 
by adjusting the magnetic field. 
The SPMHD analogue of ([13]) is 



PoPa 



dB a 

dt 



P0PaC 2 h dt 



= 0. (21) 



where the i/j energy term §T7\ has been used. Inserting ([2]), 
with no damping and (20) as the operator choice for V • B, 
produces 

EB a / dB a 
p p a \ dt /ip 

-^rn a ^^^m b (B a -B b ).\7 a W ab (h a ). (22) 

Splitting the RHS into two halves, performing a change of 
summation indices on the second half, and recombining, we 
can obtain 



dBa 

dt 



-p a ^ m b [7^2^ aW ab (h a ) 



n b p 2 b 



VaWab(h h )]. (23) 



This symmetric form for is the same as the gradient 
operator in the momentum equation. Alternatively, if this had 
been used as the operator for V B, then the differenced 
derivative operator would be imposed for V^- The occurrence 
of conjugate operators in SPH has been previously noted |T5| . 



The energy change due to damping can be written as 

'dE\ fd^ a " 



dt J 



m„ 



damp 



P0PaC% 



dt 



damp 



17l n 



PoPaC 2 h r' 



(24) 



which is negative definite. This guarantees energy may only 
be removed. 

Finally, since the additional ^(V-v) term introduced to the 
i/j evolution equation is derived using the continuity equation, 
the form for V • v should be that as in the SPH continuity 
equation. Hence, 



^a(V-V a ) 



m b (v a -v b )-V a W ab (h a ). (25) 



III. Velocity divergence cleaning for weakly 

COMPRESSIBLE SPH 

A. Weakly compressible SPH 

A common method for modelling incompressible fluid 
behaviour with SPH is to use a stiff equation of state with 
the standard Lagrangian SPH formulation. This sacrifices true 
incompressibility for simplicity of implementation. However, 
this does not imply computational efficiency as the high speed 
of sound (~ 10 x maximum fluid velocity as a minimum) 
necessitates small sized time steps for stability. Using the 
equation of state 

7 \ 



P = 



C 2 sP0 




(26) 



where po is the reference density of the fluid and c s is sound 
speed, this typically results in density variations of ~ 1% (2). 
The equations of motion which are solved are 



dt 



r = "E 



m b 



Pt 



^ ) ^oW ab . 

pI 



(27) 



In this case, we evolve the density using the SPH equivalent 
of the continuity equation, 

^ = - E m ^ v « - ^) V a W ab , (28) 

b 

rather than by summation. The smoothing length of the 
particles is held constant, calculated according to ([5]). 

B. Hyperbolic divergence cleaning for the velocity field 

Since the continuity equation relies on V • v to evolve 
density, minimising this quantity should lead to improvements 
in the representation of incompressibility. We now construct a 
formulation of divergence cleaning suitable for the velocity 
field. The cleaning equations to be solved are modified to 
become 

dv VV> 

~dt ~~ ~p~' 
d^ 2 w ^ 



(29) 



dt C ^ V ' V r 



(30) 
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As the intended application is for incompressible fluids, we 
assume throughout this section that the density is uniform and 



constant. Equations ( [29] ) and ( [30] ) still combine to produce the 
damped wave equation of ([3]). 



We follow a procedure in step with that of Sec. II-B The 
total energy of the velocity-cleaning subsystem is 



E 



I 



pdV, 



and by the constraint of energy conservation 



dE 
~dt 



F ~dt 



-X 



dtp' 
dt 



pdV = 0. 



Inserting (29 1 and (30 1 yields 



2 ^ 

-v XChPV ■ v 

P 



pdV = 0. 



Integrating the first term by parts, we obtain 

/ - - XC 2 h P (V • v)pdV + / Vv • ds = 0, 
J LP J J s 



which leads to \ = tp/c^p and hence 



i, 2 



(31) 



(32) 



(33) 



(34) 



(35) 



Hp 2 ' 

C. Discretised hyperbolic velocity divergence cleaning 

With the appropriate energy term for this cleaning system, 
the constrained SPH implementation may be constructed. We 
clean using the same V • v operator as in the continuity 
equation, that is, 



V • v a = — - ^ m b v ab • V c 

Pa , 



The SPH discretised version of ([31} is 



E = ^m a 



HA 



(36) 



(37) 



Differentiating with respect to time and using ( [30} and ( [36] ), 
we obtain 



m a v a — — 2^ — — 1^ m b v ab • V a W ah . (38) 



dt 



Pi 



By splitting the RHS into two halves, swapping summations 
on one half, then combining, it is concluded that 



dt 



(39) 



As before, conjugate operators for V • v and become 
imposed. In addition to exactly conserving energy, this form 
for also conserves momentum. 



E 

E 

x 

CO 

E^ 
CD 
> 



_ 



100000 
10000 
1000 
100 ^ 
10 

1 

0.1 
0.01 - 
0.001 
0.0001 

100000 
10000 
1000 
100 



Hyperbolic 
Hyperbolic/Parabolic 



Time 



E 

E 

x 

CD 

E 
CD 

> 



10 T- 
1 



0.1 
0.01 
0.001 
0.0001 



1 1 

Hyperbolic 
Hyperbolic/Parabolic 



4 6 
Time 



10 



Fig. 1. Maximum divergence for the free boundary test for the non- 
conservative formulation (top) and the new constrained divergence cleaning 
(bottom). For the non-conservative case, the hard boundary edge acts like 
an amplifier of divergence causing exponential growth. With the constrained 
formulation, the interaction with the boundary is treated correctly and remains 
stable. 



IV. Tests 
A. Static cleaning test: free boundaries 

The constrained cleaning methods ability to handle free 
surfaces is investigated by considering a disc of fluid with open 
boundary conditions. For this test, the full SPMHD equations 
are not solved so that the fluid retains its shape. Instead, only 
the cleaning subsystem of equations are utilised. 

The fluid is contained within a disc of radius R = 1 
composed of 1976 particles placed on a cubic lattice. The 
initial magnetic field is B z = 1 / V^tt, with a perturbation in 
the x-component of the field of the form 



B r = 



1 



(r/r ) 8 -2(r/r ) 4 



1 



r < r , (40) 
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Fig. 2. The density (top row), magnetic pressure (middle row), and divergence of B (bottom row) in the Orszag-Tang vortex at t = 1.0 comparing the 
control case (far left), including artificial resistivity (centre left), evolving the magnetic field using Euler Potentials (centre right), and applying the constrained 
divergence cleaning method (far right). 



centred on a region in the middle of the disc of radius 
ro = 1 / \/8. The density is uniformly p = 1 with zero velocity 
field. 

The maximum divergence error over time for the previ- 
ous implementation along with the new constrained imple- 
mentation of divergence cleaning is shown in Fig. [T] Both 
cases show undamped (purely hyperbolic) and damped (hy- 
perbolic/parabolic) cleaning. For the previous implementa- 
tion, once the divergence waves reach the hard edge of the 
disc, it causes divergence (and magnetic energy) to increase 
exponentially. This behaviour also occurs across jumps in 
density. However, the constrained cleaning method models the 
boundary interaction correctly. 

B. Orszag-Tang Vortex 

The constrained cleaning method is compared against arti- 
ficial resistivity and Euler Potentials using the Orszag-Tang 
vortex test problem. This problem has been widely used 
as a test of MHD codes because of its complex dynamics, 
consisting of several classes of interacting Shockwaves. 



The problem is set up in a box with dimensions x, y G [0, 1] 
with periodic boundary conditions. The initial gas state is set 
to p = 25/(36tt), P = 5/(12tt), 7 = 5/3, with velocity 
field v = [— sin(27n/), sin(27nr)]. The initial magnetic field 
is B = [— sin(27n/), sin(47nr)]. All examples presented use 
512 x 590 particles initially arranged on a hexagonal lattice. 

Results are obtained for five cases: i) no divergence control, 
ii) artificial resistivity, iii) Euler Potentials, iv) divergence 
cleaning, and v) divergence cleaning plus resistivity. Fig. [2] 
shows the density, magnetic pressure, and divergence through- 
out the system at t = 1 for the first four cases. Significant 
divergence is present in the magnetic field for the control 
case, which is reflected by small disturbances in the density 
and magnetic pressure. Artificial resistivity and Euler Poten- 
tials have an order of magnitude lower divergence error by 
comparison (Fig. [3]). Applying divergence cleaning produces 
significantly improved results. Divergence throughout the sys- 
tem is negligible, with average divergence error reduced by 
two orders of magnitude in comparison to the control case, 
down to ~ 0.1%. 
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Fig. 3. Average divergence error as a function of time in the Orszag- 
Tang vortex. Test cases included are: no divergence control, using artificial 
resistivity, employing Euler Potentials, applying divergence cleaning, and 
divergence cleaning plus resistivity. Divergence cleaning provides an order of 
magnitude reduction in divergence error over resistivity or Euler Potentials. 




Fig. 4. Column density along the y-axis of the star formation problem at 
t = l.ltff. The majority of the gas has been flattened to form an accretion 
disk about the protostar. In the left panel, magnetic divergence has grown too 
large and disrupted the system. In the right panel, divergence cleaning has 
been applied, stabilising the evolution of the system. 



Fig. 5. Average divergence error as a function of time for the star 
formation problem. Applying divergence cleaning reduces error by an order 
of magnitude. 
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C. Gravitational collapse of a magnetised molecular cloud 
core 

Our final test is drawn from our intended application: 
simulations of star formation that involve magnetic fields (16). 
These simulations follow (13), where an initial one solar 
mass sphere of gas with uniform magnetic field in the z- 
direction and in solid body rotation contracts under self-gravity 
to form a protostar with surrounding disc. However, at times 
near peak density, the magnetic field in the dense central 
region becomes strong and can produce high divergence errors. 
This has limited the range of initial magnetic field strengths 
which could be simulated, as if the divergence grows too 
large, the tensile instability correction term ^ injects enough 
momentum into the system to erroneously eject the protostar 
out of its disc 1 17 ]. Thus, this simulation proves an excellent 
demonstration of the capabilities of the constrained hyperbolic 
divergence cleaning method to reduce divergence errors in 



Fig. 6. Total linear momentum for the star formation problem. Once the 
collapse reaches peak density (t ~ 1), a sharp increase in momentum 
occurs due to divergence errors. The divergence cleaned system reduces this 
momentum spike by two orders of magnitude. 

realistic, 3D simulations. 

The sphere of gas has radius R = 4 x 10 16 cm with uniform 
density p = 7 A3 x 10 -18 g cm -3 . A barotropic equation 
of state is used, as described in (13). The magnetic field 
strength is set to give a mass-to-magnetic flux ratio of 5 
times the critical value for magnetic fields to provide support 
against gravitational collapse. To avoid edge effects with the 
magnetic field, the sphere is embedded in a periodic box of 
length AR containing material surrounding the sphere set in 
pressure equilibrium with density ratio 1:30. This test uses 
only a minimal amount of resistivity, with a# G [0, 0.1]. Self- 
gravity is simulated using a hierarchical partitioning tree, with 
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Fig. 7. Snapshots of the oscillating water drop test. The circular drop has an initial velocity which squeezes it into an elliptical shape along the y-axis. A 
radial force is present which halts the expansion of the drop, then contracts it to its original shape before expanding along the opposite axis. This behaviour 
repeats causing the drop to oscillate alternately along the two axes. 



gravitational force softening using the SPH kernel as described 
by p8) . The free fall time is ~ 24000 years. A sink particle is 
inserted once the gas density surpasses p S i n k = 10 -10 g cm -3 , 
and accretes particles within a radius of 6.7 AU. 

Fig. [4] shows column density comparisons of simulations 
with (right) and without (left) divergence cleaning at t = 1.1 
free fall time, showing that drastic improvements to the 
results are obtained by incorporating divergence cleaning. The 
protostar remains stable in its disc and a helical shaped jet is 
launched from the centre fl6| . The average divergence error 
is reduced by an order of magnitude (Fig. [5]), and this leads to 
a corresponding improvement in the momentum conservation 
of roughly two orders of magnitude (Fig. [6]). 

D. Oscillating water drop test 

To investigate the effectiveness of our velocity cleaning 
algorithm, it is applied to an oscillating elliptic water drop. 
The water drop is initially circular and is free standing. A 
radial force is exerted upon it, and with an initial velocity 
which is compressional along one axis, the drop oscillates, 
squeezing alternately along the two axes. This behaviour is 
demonstrated in Fig. [7] 

The drop is modelled using the weakly compressible ap- 
proximation (equations ( [27] ) and ( [28] ) with (26) as the equation 
of state). The reference density is po = 1000 kg m -2 , and the 
initial velocity field is v = [— 100x, 100?/]. The radial force 
is — 100 2 r. The drop has radius R = 1, and a total of 1976 
particles are used arranged on a square lattice. 

The evolution of the drop is tracked until t = 0.1, ap- 
proximately two oscillation periods. Fig. [8] shows the average 
velocity divergence of the system as a function of time for both 
the cleaned and uncleaned systems. Applying cleaning reduces 
the average divergence by nearly an order of magnitude, 
similar to results obtained for magnetic field cleaning. This 
leads to a reduction in maximum density error by a factor of 
two (Fig. [9]). The dissipation of kinetic energy by the cleaning 



100 



algorithm is insignificant, as shown in Fig. 10 



V. Conclusion 

In this paper, SPH formulations of Dedner et al's hy- 
perbolic/parabolic divergence cleaning for the magnetic and 




0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 
Time 

Fig. 8. Average V • v of the elliptic water drop test. Average velocity diver- 
gence is reduced by approximately an order of magnitude when divergence 
cleaning is applied. 




Control 

Divergence Cleaning 




V V V V' 



0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 
Time 

Fig. 9. Maximum density variation during the elliptic water drop test. 
Applying divergence cleaning to the velocity field reduces density changes 
from the reference density by ~ 0.5. 
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Fig. 10. Total kinetic energy of the elliptic water drop test. No significant 
discrepancies exist between the control and divergence cleaned tests. 



velocity fields have been presented. The algorithm is attractive 
because it is computational inexpensive and easy to implement 
in existing codes. For SPMHD simulations in particular, it rep- 
resents a path forward for maintaining the magnetic divergence 
constraint without the drawbacks associated with using Euler 
Potentials or artificial resistivity. 

Our method was derived by considering the energy con- 
tained in the ip field as part of the system total energy. With 
this contribution included, it is possible to construct SPH 
implementations which conserve energy and, in the case of 
the velocity field, momentum. This stabilises the algorithm 
across density jumps and at free surface boundaries. Results 
obtained find an order of magnitude reduction in average 
level of divergence for all tests of magnetic and velocity field 
cleaning. For 3D astrophysical star formation problems, this 
reduction in magnetic divergence has improved momentum 
conservation by two orders of magnitude. When the velocity 
field is cleaned in weakly compressible SPH simulations of an 
oscillating water drop, density errors are reduced by half with 
negligible kinetic energy dissipation. 

Given the performance of the algorithm on complicated as- 
trophysical applications, the door to study other astrophysical 
problems is now open. Though the results of velocity cleaned 
weakly compressible SPH are encouraging, additional work is 
required, in particular, for cases involving boundary particles. 
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